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LIST  OF  SYMBOLS  AND  ABBREVIATIONS 

A  Arbitrary  constant  in  the  stream  function 

a  Earth's  radius 

F  Zonal  flux  term 

g  Flux  term  in  ^-equation  of  motion 

g  Flux  term  in  ri-equation  of  motion 

G  Meridional  flux  term 

I  Number  of  points  around  latitude  circle 

i  Grid  index  in  the  ^-direction 

j  Grid  index  in  the  rpdirection 

k  Vertical  grid  index 

m  1/ (a  cos  <{>) 

mb  Millibars 

NACA  National  Advisory  Committee  on  Aeronautics 

n  1/a 

P  Pole  on  the  ij  index  system 

• 

S  Area-pressure  weighted  vertical  velocity 

t  Time 

u  Zonal  wind 

v  Meridional  wind 

n  Meridional  coordinate  of  the  curvilinear  coordinate  system 

An  Distance  increment  in  the  meridional  direction 

X  Longitude 

^  Zonal  coordinate  of  the  curvilinear  coordinate  system 

A^  Distance  increment  in  zonal  direction 


rr  Terrain  pressure 

\\  Area-weighted  terrain  pressure 

a  Dimensionless  vertical  coordinate 

Ac  Vertical  increment  in  the  sigma  coordinate  system 

a  Measure  of  vertical  velocity 

♦  Latitude 

Q  Angular  velocity  of  the  earth 


I.   INTRODUCTION 

Most  global  primitive  equation  models  now  in  use  employ  spherical 
coordinates  (see  Haltiner  and  Williams  [1975]).   Flow  crossing  the  pole 
presents  a  potential  problem  in  this  coordinate  system.  Mihok  and  Kaitala 
(1976)  have  described  a  global  prediction  model  which  is  in  the  last  stages 
of  development  at  the  Fleet  Numerical  Prediction  Central.   McCollough  (1974) 
in  testing  an  earlier  version  of  this  model  observed  the  development  of  large 
gradients  in  the  surface  pressure  near  the  pole  when  a  particular  real  data 
set  was  used.   Maher  (1974)  examined  this  problem  more  closely  by  using  an- 
alytic initial  data  which  gave  a  strong  f low  across  the  pole.   His  solutions 
showed  that  the  pressure  field  was  spuriously   disturbed  by  the  finite  dif- 
ferencing near  the  pole.  He  was  able  to  reduce  the  effect  by  various  types 
of  smoothing. 

In  this  report  we  will  test  a  somewhat  different  finite  difference 
scheme  and  we  will  test  a  procedure  for  controlling  the  problem.   The  Navy 
Environmental  Prediction  Research  Facility  is  now  testing  the  global  pre- 
diction model  which  was  developed  by  Arakawa  and  Mintz  (1974).   In  this 
study  we  will  employ  the  model  which  was  described  by  Monaco  and  Williams 
(1975).   This  model  is  a  slightly  simplified  version  of  the  Arakawa  model 
and  it  differs  from  the  FNWC  model  mainly  with  respect  to  the  spatial 
staggering  of  the  variables. 


II.   FINITE  DIFFERENCE  EQUATIONS 

The  basic  differential  and  difference  equations  used  in  this  study 
are  given  in  Monaco  and  Williams  (1975).   Here  we  will  reproduce  only 
those  difference  equations  which  have  been  changed.   Fig.  1  shows  the 
arrangement  of  variables  in  the  model.   The  list  of  symbols  gives  the 
various  definitions.   The  finite  difference  continuity  equation  has 
the  following  form: 


dp*- 1  i  +  Fk       vk  ,  rk       rk 

Ot         1+2, J      l-i,J      l,J+2      l,J-2 


+  1T(Sk+1  -  Sk"b  =  0.  (2.1) 

Aak  X'J    X'J 


The  mass  fluxes  are  defined 


'i+i  -  =  4<«^>  >4i  .(tt.     +  TT.  .),  (2.2) 

1+2  >  J         n  l+t,J   1+1  >  J      IjJ 


x,j+2-        m  l,j+2   1,J+1      1,1 


where 

I*. 


j  ~  TTi,j(i^r) .  . 
i»j 


i»j   »!i,j  i,j 


The  bar  in  (2.2)  represents  zonal  smoothing  as  described  inMW. 


Hereafter  we  will  refer  to  this  reference  as  MW. 


n    CENTERED 
NOTAT ION 


u  CENTERED 
NOTAT ION 


Figure  1.   An  example  of  the  notation  used  to  describe  the  finite  dif- 
ference equations.   The  continuity  equation  is  described  on  a  "ir-centered" 
grid  in  which  the  F  and  G  symbols  are  flux  calculations  in  their  re- 
spective directions.   The  "u-centered"  grid  is  an  example  used  to  de- 
scribe the  ^-component  of  the  equation  of  motion  where  F      and  g  are 

also  flux  calculations. 
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The  left  hand  side  of  the  zonal  momentum  equation  is: 


at  Ui,j  i,j/   «*  i+|,j  i+l, j    x,j 


-  F.  i  .(u    +  u     )   +  g   .ji(u   .,-  +  u   )   -  g    i(u    +  u     )  J 

1~2,J        1>J  1-1»J  l,J+2        IjJ+1  1,J  1»J~2        1»J  1,J-1 


l     ir*u,k+l.k+2         k      s      •u.k-l,  k        ,      k-2x-,  ..    .  , 

h r  t  S.    .      (u,    .  +  u.    . ;-  S.    .      (u.    .  +  u.     JJ  (2.4) 

k  *     I,j         i,j         i,j         i,j         i,j         i,j 


where 

fli4'ki«  +iii-i,j+i  +  M,j-i  +m4,j-i  +2(M,j  +I\i4,3)]'  (2-5) 

^i,j  "  K^.J+l  +  Sl4;j+1  +  8l4i>1  +  S^.,  +  2(S^;.  +  Sl4>.)].    (2.6) 


The  above  definitions  (2.5)  and  (2.6)  are  more  complicated  than  those  in  MW 
and  they  are  the  same  as  those  used  by  Arakawa  and  Mintz  (1974) .  The  quan- 
tities F  and  e  are  defined  as  follows: 


««*.j s  y<f-4.*i +  °t4.i  +  FU,j.i>>  <2-7> 


S?  5JA  £  7(G*aJL  •  +  G*^k  -u.1  +  G*  i  .  +  G*  i  ..,),  (2.8) 

l.J+S"    *    1+2, J      1+2, J  +  l      3.-5,  J      1-2,  J  +  l 


where 


Ft   =i(F  ^   +  F  ,   ),  (2.9) 

X,J         1+2,  J      1**2,J 


G*  .  =  |(G.  .  ,  +  G.  .  x).  (2.10) 

i,J      1,3+1    i,J-f 


The  left  hand  side  of  the  meridional  momentum  equation  takes  the 
same  form  as  (2.4)  when  u  is  replaced  by  v.  The  other  quantities  are 
defined  as  follows: 


fY^i  •  =   f(F*_ .  .j,  +  F*   .  +  F*  .  i  +  F*    .  i),  (2.11) 

i+f, J   4   i+l, j+t    i,J+£    i,j-§    1+1, j-t 


5V    -J.  -  fCG*.,     .j,  +  2G*    .   ,    +  G*   .     .    .),  (2.12) 

i»J+t        4      i+l, J+5  i,J+t  i-lfJ+t 


v  1 


TFi.j  s  sfi+ij+i  +jfi-i,j+i  +/Ti-i,j-i  +Ti+i5j-l  +  2(ili,j+l  +'|U,j4)]'    (2a3) 

•  1*  •  •  •  •  * 

SY  .  =   ^[S.,,  ._,!  +S.  .  .,i  +  S.  .  .  i  +S.^  .  !  +  2(S.  .i  +  S.  .  i)].    (2.14) 
i,j    8  l+l, j+t   l-l, j+l    i-l, J-f   i+l, J-t      i,J+f    i,J-t 


The  definitions  (2.13)  and  (2.14)  are  more  complicated  than  those  in  MW 
and  they  are  the  same  as  those  used  by  Arakawa. 

Arakawa  and  Mintz  (1974)  derived  the  expressions  (2.5)  and  (2.6)  for 
(I  and  S  at  the  u  points  in  the  following  way.  The  u  field  in  (2.4)  is  set 
equal  to  a  constant.  Then  the  continuity  equation  (2.1)  is  averaged  from 
surrounding  points  in  such  a  way  that  the  average  has  the  same  form  as 
(2.4)  with  u  =  constant.  This  is  a  reasonable  requirement  and  it  is  nec- 
essary for  energy  conservation.  This  average  gives  the  expressions  (2.5) 
and  (2.6).  The  definitions  (2.13)  and  (2.14)  are  derived  in  the  same  way 
by  setting  v  =  constant  in  the  appropriate  equations. 

These  difference  expressions  must  be  modified  near  the  pole  because 
some  of  the  quantities  are  not  defined  at  the  pole.   Fig.  2  shows  the 
points  near  the  pole  which  are  used  to  evaluate  the  pressure  change  at 
the  pole.   The  pole  is  composed  of  I  points  which  are  denoted  by  (i,  P) . 
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Figure  2.   Each  index  i,P  in  the  continuity  equation  is 

represented  by  the  shaded  area,   tt  at  the  poles 
can  only  change  as  a  result  of  G.   The  thermo- 
dynamic equation  and  vertical  velocity  are  treated 
in  the  same  manner. 


The  continuity  equation  applied  at  each  of  these  points  is  written: 


V  «*.***  £33  -  £i>  -  «■ 


After  each  time  step  the  surface  pressure  at  the  pole  (P)  is  set  equal 
to  the  average: 

'-liji.p  .  (2-l6) 


If  we  sum  (2.18)  in  the  vertical  and  use  (2.16)  we  can  see  that  the  pres- 
sure at  the  pole  changes  in  proportion  to  the  net  mass  flux  across  the 
latitude  circle  at  j  =  P-§-. 

Fig.  3  shows  the  quantities  which  are  needed  to  predict  u  at  j  =  P-l. 
The  left  hand  side  of  the  zonal  momentum  equation  at  j  =  P-l,  takes  the  form: 

lt^^iXi+*^,M(ui>i  +ui+i.p-i)k 

-Fi'4,P-l(Ui,P-l  +  Vl,P-l)k  "  gi,P-3/2(ui,P-l  +  Ui,P-2)k] 

,  1   lr*u,k+l  k+2      k        lu,k-l.  k         k-2   ,       ,   7 

+  T'[S   K,p-i  +  ui,P-i  >  "  s    (ui,P-i  +  ui,P-2)] '     (2'17) 

where 

Ff4,p-i  n(WH,p-i  +  rt-i,p-2»  '  (2a8) 
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Figure  3.   The  polar  modification  of  the  u  equation  of  motion 
F        and  g   are  flux  terms.   The  shaded  portion 
represents  the  area  associated  with  each  variable 
u. 
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¥ i.p-i s  iiP  +  iC-§,p-i  +/Ti^,p-i)  +  l^i-i.p-2  +TTi+i)P-2)  ■   (2-19: 


*i,P-l  !  S  +  l«l-i,P-l  +  "i-A.P-l*  +  i(Vi.P-2  +  =l+i,P-2)-    (2'20) 


The  other  quantities  have  the  same  definitions  as  above.   The  relations 
(2.19)  and  (2.20)  are  different  from  those  used  by  MW  and  they  correspond 
to  the  relations  derived  by  Arakawa. 

Fig.  4  shows  the  quantities  which  are  needed  to  predict  v  at  j  =  P-|-. 
The  left  hand  side  of  the  meridional  momentum  equation  at  j  =  P-^j  takes  the  form: 


4-(nVvk)     i  +  ir[F*Vi    i(v.      1  +v    i) 
dt"    V  ;i,P4   2L  i+i,P-ir   i+l,P4    i,P"| 


■f^,p4(vi,p4  +  vi-i,p4)k  "  gI,p-i(vi,p4  +  vi,p-3/2>ki 

+  2_^,^2^+v^    _   s^-1^)   +v^4)],  (2.21) 


where 


*&,*-* -4<pt-i.p-i  +  FW'  (2-22) 


Hl.p-i  SV  s^i-i.p-i  +  TTi+i,p-i  +  2?i,P-r 


(2.23) 


sI.p-i  !  ^ +  i»i-i,p-i +  si+i.p-i  +2=i,p-1)  •  (2-24) 
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Figure  4e   The  polar  modification  of  the  v  equation  of 

motion.  F  and  g  are  flux  terms.  The  shaded 
portion  represents  the  area  associated  with  each 
variable  v. 
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The  relations  (2.23)  and  (2. 24)  are  different  from  those  used  by  MW  and  they 
are  equal  to  the  relations  derived  by  Arakawa. 

The  expressions  for  Jy    and  J[     near  the  pole  were  derived  by  Arakawa 
and  Mintz  (1974)  in  the  same  manner  as  for  other  latitudes.   For  example, 
u  is  set  equal  to  a  constant  in  Eq.(2.17)  and  the  continuity  equation  is 
averaged  to  achieve  the  same  form.   However,  the  derivation  used  (2.15), 
but  does  not  use  the  average  condition  (2.16). 

III.   NUMERICAL  SOLUTION  WITH  FLOW  OVER  POLE 

In  this  section  we  will  examine  a  numerical  solution  which  was  ob- 
tained with  the  difference  equation  described  in  section  2.   The  initial 
conditions,  which  are  similar  to  those  used  by  Holloway,  Spelman  and 
Manabe  (1973),  contain  strong  flow  across  the  poles.   The  initial  con- 
ditions are  taken  from  the  barotropic  Rossby  wave  solution  developed  by 
Haurwitz  (1940) .   These  initial  data  lead  to  an  exact  solution  for  baro- 
tropic horizontal  motion  (Neamtan  [1946]).   However,  in  our  model  the  di- 
vergence will  not  remain  zero,  but  the  exact  solution  should  be  close  to 
the  Haurwitz  solution. 

The  initial  zonal  wind  is  given  by 


A         2       2 
u(\,4>,o,0)  = sin\[cos  *  -  sin  4>]  ,  (3.1) 

EL 


and   the  meridional   component   is   given  by 


v(\,*,a,0)   =  -  cosX    sin*.  (3.2) 

ci 
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These  fields  describe  wavenumber  1  planetary  flow  with  no  mean  current, 
The  initial  surface  pressure  is  given  by 


rr(X,*,t)  =  tt  „  +  K[(^)2(2cos4*  -  cos2*  -2) 
ave      /a 


+  9^cos  $(5-4  cos2*)  sinX  +  (~)2  cos2*  (3  -  2cos2*)cos  X  ]  .    (3.3) 
3  2a 


Here  tt    is  the  mean  surface  pressure  and  K  is  a  constant  of  proportion- 
ave 

ality.   This  solution  for  the  pressure  field  was  obtained  by  Phillips  (1959) 
and  is  equivalent  to  the  solution  of  the  nonlinear  balance  equation.   The 
initial  temperature  field  is  a  function  of  height  only,  and  it  is  given  by 
the  NACA  standard  atmosphere. 

Fig.  5  shows  the  initial  surface  pressure  field  in  the  vicinity  of  the 
pole.  Note  the  strong  geostrophic  flow  directly  across  the  pole. 

The  tests  described  in  this  report  were  carried  out  with  a  two-level 
model,  with  45  points  between  the  poles  and  10  points  in  the  east-west 
direction.  A  time  step  of  6  minutes  was  used.   No  heating,  friction  or 
surface  topography  were  included. 

Fig.  6  shows  the  surface  pressure  prediction  of  t  =  54  hours  which  was 
made  from  the  above  initial  conditions.   The  basic  pattern  has  rotated  con- 
siderably from  its  initial  orientation  because  the  wave  number  1  initial 
state  excites  a  Rossby  wave  which  moves  rapidly  westward.   However,  in  the 
vicinity  of  the  pole  the  isobars  have  become  quite  distorted.   In  the  non- 
divergent  solution  obtained  by  Neamtan  (1946)  the  initial  disturbance  rotates 
at  constant  angular  velocity  without  distortion.  Although  in  our  model  the 

15 


Figure  5.   The  initial  surface  pressure  field  in  the  vicinity  of  the  pole. 
The  location  of  the  pole  is  indicated  by  the  small  circle  in  the  center. 
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Figure  6.   The  surface  pressure  field  in  the  vicinity  of  the  pole  at 
t  =  54  hours. 
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divergence  is  not  zero,  it  is  expected  that  the  general  behavior  should 
be  similar  to  the  behavior  of  the  nondivergent  solutions.   Thus  the  dis- 
tortions which  appear  in  Fig.  6  near  the  pole  are  of  a  computational 
nature.   In  fact,  the  numerical  solution  "blew-up"  at  about  t  =  75  hours. 

In  order  to  consider  this  problem  further,  let  us  examine  the  V 
component  of  the  velocity  near  the  pole.   In  the  upper  part  of  Fig.  7  is 
the  initial  v  component  at  <j>  =  88  ;  the  v  component  at  $  =  84   is  almost 
identical.   In  the  lower  portion  of  Fig.  7  are  the  v  fields  at  the 
2  latitudes  for  t  =  54  hours.   The  fields  are  somewhat  out  of  phase  with 
each  other  and  the  component  nearest  to  the  pole  is  greatly  amplified. 
After  t  =  54  hours  the  maximum  v  component  continues  to  grow  until  the  num- 
erical solution  is  completely  destroyed. 


IV.   CONTROL  OF  POIAR  PROBLEM 

In  Fig.  7  it  was  seen  that  the  v  field  along  the  row  of  points  next  to 
the  pole  becomes  very  large  and  irregular  as  the  numerical  solution  becomes 
unrealistic.   It  appears  that  the  geostrophic  adjustment  process  cannot 
operate  effectively  along  the  row  of  points  next  to  the  pole  ($  =  88   in 
our  experiments).   For  example,  if  the  v  component  is  too  large  at  one 
point,  then  it  would  modify  the  polar  value  of  the  pressure,  which  would 
change  the  pressure  gradient  at  the  original  points.   However,  this  pro- 
cess cannot  operate  freely  because  the  polar  pressure  changes  in  propor- 
tion to  all  the  points  of  $  =  88  ,  not  just  one  (see  Eqs.  (2.15)  and  (2.16)) 
It  is  true  that  various  quantities  near  the  pole  such  as  (2.22)  and  (2.23) 
and  (2.24)  were  derived  by  assuming  consistency  between  the  momentum  change 
at  4>  =  88  and  the  pressure  changes  at  surrounding  points.   But  this  proof 
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36CT 


Figure  7.   The  v  field  as  a  function  of  longitude  at  t  =  0  and  t  =  54  hours 
at  the  latitudes  indicated. 
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is  not  strictly  valid  because  it  assumes  different  values  of  tt-  p 
(see  Eq.  (2.15))  for  different  longitudes  whereas  (2.16)  is  enforced  at 
every  time  step. 

In  order  to  avoid  this  problem  it  seems  appropriate  to  keep  only  that 
portion  of  the  v  field  at  4>  =  88   that  can  achieve  geostrophic  balance.   This 
restricts  the  v  field  to  just  wave  numbers  0  and  1.  With  wave  number  0  we 
have  the  proper  interaction  with  the  polar  pressure,  since  v  is  independent 
of  longitude.   The  velocity  vector  is  the  same  on  both  sides  of  the  pole  with 
wavenumber  1,  so  that  the  polar  pressure  should  not  change.   Thus  after  every 
time  step  we  set 


1         2    1 


i,p4"if4V4     n  =  i  vi,p4     v  i  ;  '       v  i 
2uf  v     „       _e-,w£ni*    uin^M 


+  i(iSi  vi,p-|sin(i^    )sin(  i  >•  (4-1} 

Fig.  8  shows  the  surface  pressure  field  at  t  =  54  hours  when  (4.1)  is 
used  every  time  step.   In  this  case  the  isobars  are  still  distorted  near 
the  pole,  but  not  nearly  as  much  as  in  Fig.  6.   In  Fig.  9  we  see  the  pre- 
dicted v  fields  at  the  2  highest  latitudes  for  t  =  54  hours.   These 
fields  are  much  closer  to  the  initial  fields  than  those  shown  in  Fig.  7. 
There  is  a  small  increase  in  the  velocity  at  $  =  88  and  both  fields  are 
quite  smooth.   Of  course  v.   i  is  forced  to  be  smooth  by  (4.1).   Fig.  10 
shows  the  surface  pressure  field  at  t  =  108  hours.   It  is  actually  smoother 
than  at  t  =  54  hours,  although  the  later  solution  does  have  some  tilt.   It 
should  be  noted  that  the  total  disturbances  amplitude  has  decreased  by  about 
20%  at  t  =  108  hours.   This  is  apparently  due  to  the  Euler -backward  restart 
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Figure  8.   The  surface  pressure  field  at  t  =  54  hours  obtained  with  the 
modified  v  structure. 
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Figure  9.   The  v  field  as  a  function  of  longitude  obtained  with  the 
modified  v  structure  for  t  =  54  hours  and  t  =  108  hours » 
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Figure  10.   The  surface  pressure  field  at  t  =  108  hours  obtained  with 
the  modified  v  structure. 
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which  is  used  in  the  model  every  30  minutes  and  to  the  high  frequency  of 
the  Rossby  wave  in  this  case.   The  v  field  near  the  pole  is  smooth  and 
has  a  small  amplitude  as  may  be  seen  in  the  lower  portion  of  Fig.  9. 


V.   CONCLUSIONS 

It  has  been  shown  that  the  global  model  developed  by  Monaco  and 
Williams  (1975)  is  inaccurate  when  there  is  flow  over  the  pole,  and  in 
fact  the  model  may  become  unstable.   This  instability  was  attributed  to 
the  difficulty  in  achieving  geostrophic  adjustment  near  the  pole.   This 
is  because  the  polar  pressure  is  changed  by  the  average  mass  flux  across 
the  row  of  points  nearest  to  the  pole.   Thus  the  polar  pressure  reacts 
only  weakly  to  the  v  at  a  single  point  adjacent  to  the  pole. 

The  instability  was  eliminated  when  the  v.    i  field  was  replaced 
by  its  average  plus  the  first  wave  in  longitude.   The  resulting  integra- 
tions were  stable  and  the  fields  were  smooth,  although  some  distortion  still 
occurred  near  the  pole.   Sadourny  (1975)  has  formulated  a  global  model  in 
cylindrical  coordinates  which  conserves  mean  square  potential  vorticity. 
His  solutions  with  wave  number  1  appear  to  be  quite  good,  although  he  does 
not  show  the  detail  near  the  pole. 

The  technique  developed  in  this  report  may  be  useful  for  controlling 
the  polar  problem  in  other  global  models.   It  could  be  used  in  the  model 
described  by  Mihok  and  Kaitala  (1976)  which  is  now  under  development  at 
the  Fleet  Numerical  Weather  Central.   In  this  application  the  fields  u,  v, 
n  and  T  would  all  have  to  be  treated  since  they  all  occur  on  the  same 
latitude  circle. 
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